A universal variational quantum eigensolver for non-Hermitian systems

Many quantum algorithms are developed to evaluate eigenvalues for Hermitian matrices. However, few practical approach exists for the eigenanalysis of non-Hermintian ones, such as arising from modern power systems. The main difficulty lies in the fact that, as the eigenvector matrix of a general matrix can be non-unitary, solving a general eigenvalue problem is inherently incompatible with existing unitary-gate-based quantum methods. To fill this gap, this paper introduces a Variational Quantum Universal Eigensolver (VQUE), which is deployable on noisy intermediate scale quantum computers. Our new contributions include: (1) The first universal variational quantum algorithm capable of evaluating the eigenvalues of non-Hermitian matrices—Inspired by Schur’s triangularization theory, VQUE unitarizes the eigenvalue problem to a procedure of searching unitary transformation matrices via quantum devices; (2) A Quantum Process Snapshot technique is devised to make VQUE maintain the potential quantum advantage inherited from the original variational quantum eigensolver—With additional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(log_{2}{N})$$\end{document}O(log2N) quantum gates, this method efficiently identifies whether a unitary operator is triangular with respect to a given basis; (3) Successful deployment and validation of VQUE on a real noisy quantum computer, which demonstrates the algorithm’s feasibility. We also undertake a comprehensive parametric study to validate VQUE’s scalability, generality, and performance in realistic applications.

Eigenvalue calculation forms the backbone of a multitude of scientific and engineering disciplines, many of which necessitate the evaluation of eigenvalues for non-Hermitian matrices.Particularly, in power systems' dynamic analysis, the eigenvalues of non-Hermitian matrices are critical for determining the system's small signal stability and oscillation modes [1][2][3] .As the grid-connected renewable energy resources swiftly increase and reach a total capacity of over 3.4 TeraWatts globally 4 , most of which being wind farms and photovoltaic solar arrays, these intermittent, inertia-less renewable electricity injections, in addition to the weather and climate impacts as well as unforeseeable faults, are exposing unprecedented challenges on the stability and resilience of power grids worldwide.Consequently, there is an urgent need to evaluate eigenvalues for extremely large non-Hermitian matrices arising from power systems so as to enable fast assessment and assurance of power system stability and resilience.
The QR algorithm, underpinned by Schur's decomposition, is the most commonly employed method for calculating eigenvalues 5 .However, it demands a computational complexity of O(N 3 ) , posing a significant chal- lenge for analyzing large networks due to its scalability limitations.In pursuit of greater scalability, extensive efforts have been devoted to the development of iterative Krylov subspace eigenvalue algorithms, such as the Arnoldi algorithm 6 .By leveraging the sparsity property, these methods can efficiently compute a limited subset of eigenvalues, making them particularly the method of choice for large network analysis.
The fast development of quantum computation technologies carries a ray of light to forward this field further.For the eigenvalue evaluation problem of Hermitian matrice, the quantum phase estimation (QPE) method, for example, can find eigenvalues exponentially faster than the classical methods 7,8 ; however, the probability distributions of eigenvalues it depends on the respective overlap of the initial states with eigenvectors.Recently, by generalizing the classical power method and employing techniques in the famous quantum linear system (known as HHL 9 ), Nghiem and Wei 10 introduce a quantum algorithm to estimate the largest eigenvalues of a Hermitian matrix, which also inherits the advantage of the HHL algorithm.Significant efforts have been expended in the last decade to extend the QPE method to the non-Hermitian matrix [11][12][13] .As shown in a recent research 14 , QPE is extended to evaluate the eigenvalues of a non-Hermitian but diagonalizable matrix, which either only has real eigenvalues or the corresponding input state could be prepared naturally.Furthermore, these methods all require www.nature.com/scientificreports/an excessive long qubit decoherence time, which is unachievable on today's noisy intermediate-scale quantum (NISQ) devices.The aforementioned problems are partially overcome by the variational quantum eigensolver (VQE) method proposed by Peruzzo et al. 15 .By enabling the quantum computer to work in concert with classical computer resources (e.g.optimizer), the long quantum process of finding eigenvalues can be broken up into a combination of multiple parameterized quantum programs (usually called ansatz).These ansatzes are shallow quantum circuits that only require a manageable coherent time.A classical optimizer assembles them globally to accomplish the eigenvalue evaluation.Such hybrid quantum-classical strategies are known as variational quantum algorithm (VQA) and have emerged as the leading technique to obtain a quantum advantage in the current and near future of NISQ era 16 .Such variational approaches work particularly well for obtaining the largest or smallest eigenvalues of a Hermitian matrix.Tremendous efforts are still being made to improve the VQE algorithms and ascertain their advantage.For instance, a large class of quantum algorithms are developed to enable VQE to search for excited states [17][18][19][20][21] .By exploiting the relation of diagonalization and majorization, Cerezo et al. 22 introduces a more amenable method for estimating the largest eigenvalues.To improve the convergence and evaluation accuracy, Harwood et al. 23 introduces a hybrid homotopy continuation algorithm.A more detailed VQE method evolution history can be found in the review papers by Cerezo et al. and Tilly et al. 16,24 .The eigenvalue algorithm is also extended to find the singular value of an arbitrary matrix 25,26 , as the singular value decomposition is realized by two unitary transformation matrices.Very recently, a variational method has been proposed 27 to evaluate a specific eigenvalue of a non-Hermitian system (This paper appeared around the same time as our manuscript was submitted.We acknowledge a reviewer for bringing it to our attention.).This method requires simultaneous optimization of both the eigenvalue and the ansatz to yield an eigenvalue estimation near a pre-specified value.However, to estimate each unique eigenvalue, the entire optimization must be performed anew, rendering the process inefficient for determining multiple eigenvalues of a large matrix.
However, the gap remains between the developed quantum eigenvalue algorithms and real-world power system (and most other engineering systems) problems as all power systems, no matter direct current (DC), alternating current (AC) or hybrid, have non-Hermitian systems matrices.The primary challenge is that the eigenvectors of a non-Hermitian matrix can be non-unitary, leading to instances where the use of unitary gatebased quantum methods fails to identify the eigenvector matrix required for matrix diagonalization.Further, power system small signal analysis requires identifying all the complex eigenvalues to determine the oscillation modes of the power network.Therefore, few deployable quantum algorithm has been established thus far to perform eigenanalysis for an arbitrary matrix.
It is also worth noting that, by using the method presented in the HHL paper 9 , a non-Hermitian matrix can be converted to a Hermitian matrix for eigenvalue evaluation.However, the eigenvalues of this Hermitian matrix are related to the non-Hermitian matrix's singular values.Since generally no relationship exists between eigenvalues and singular values, the eigenvalues of the non-Hermitian matrix are still unknown and cannot be obtained by methods designed for the latter.
It is also worth noting that a non-Hermitian matrix can be converted to a Hermitian matrix for enabling quantum computation, as introduced in the classic HHL work 9 .However, the eigenvalues of this resultant Hermitian matrix are related to the non-Hermitian matrix's singular values.Since there is generally no relationship between eigenvalues and singular values, the eigenvalues of the non-Hermitian matrix are still unknown and cannot be obtained by methods designed for the latter.
In addition, classical VQE methods only give estimated eigenvalues, while no accuracy information is directly available.However, it is desirable to quantify the estimation errors for practical applications.
To address the aforementioned challenges, we devise a novel variational quantum universal eigensolver (VQUE) which empowers eigenanalysis of a generic matrix.Specifically, our contributions are threefold: • By reexamining Schur's theory, we unitarize the eigenvalue problem into a process of searching a unitary transformation matrix.This forms a cornerstone of a practical variational quantum algorithm capable of finding eigenvalues for non-Hermitian systems.• We devise a quantum process snapshot method to determine if an N-dimensional matrix is triangular.
Compared with the brute-force approach, it reduces the computational complexity from O(N 2 ) to O(log 2 N) gates.Thus the quantum advantage of VQUE is maintained.We also propose a novel approach to implement the non-unitary operator directly within the quantum process snapshot, thereby enabling VQUE to examine any arbitrary matrix.• We define our cost function as the distance between the transformed matrix and a triangular matrix.This can be conveniently determined using the newly proposed quantum process snapshot technique.In contrast to the traditional VQE method, this cost function also doubles as a gauge for accuracy.Therefore, it enables the use of a shallower tentative ansatz in a trial-and-error process, which substantially enhances VQUE's performance in real-world applications.
VQUE's scalability, universality, and efficacy are validated on noise-free simulators and NISQ devices.

Preliminaries Variational quantum eigensolver
To establish a link between the Variational Quantum Eigensolver (VQE) algorithm and the mathematical process of identifying the smallest eigenvalue of a Hermitian matrix, this subsection succinctly revisits the VQE method.We will employ matrix and vector notation for this recapitulation, as demonstrated below.
For a Hermitian matrix H (e.g., hamiltonian matrix), its Rayleigh quotient is defined as : where the superscript † denotes complex conjugate and, without loss of generality, and when x is taken as a vector with unit magnitude then the denominator can be omitted.
It can be shown that the Rayleigh quotient reaches its minimum value R(H, � x) min = min , when x is the eigenvector corresponding to the eigenvalue min 5,28 .Based on the mathematical theory above, the VQE method decomposes the Hermitian matrix as a linear combination of tensor products of Pauli operators 15,29 : where σ s denotes the s-th unitary matrix which is formed by n k=1 σ s,k and c s = 1 2 n Tr(σ s H) is the s-th coef- ficient correspondingly; σ s,k ∈ {σ I = I, σ x , σ y , σ z } is the single-qubit Pauli operator in s-th unitary matrix for the k-th qubit.
Therefore, by applying a parameterized circuit U (θ) to a fixed initial state such as |0 • • • 0� , i.e., with the output vector � x = U (θ)� x 0 , the algorithm obtains an estimate � The estimate is iteratively optimized by a classical controller, which changes the parameter θ to minimize the expectation value x † H x (or �x|H|x� in the standard Dirac notation) obtained by measurements, which equals to: The quantum advantage is then obtained from evaluating each x 0 † U (θ) † σ s U (θ) x 0 15 on a quantum computer.In a classical computer, calculating such matrix multiplication with N dimension takes O(N 2 ) computation costs.As a comparison, a quantum computer has the capability to evaluate such terms (also known as evaluating the expectation value) by joint local measurement of each qubit, which is exponentially faster than the classical computer 8,30 .
However, all steps in the above procedure are based on the pre-assumption that the matrix H is Hermitian, whose eigenvector matrix is unitary.

Mathematical foundation
The critical aspect of enabling variational quantum algorithms to calculate the eigenvalues of a general matrix involves identifying a unitary transformation matrix that can distinctly expose the eigenvalues in a similar manner to the eigenvector matrix.In our devised variational quantum universal eigensolver (VQUE), we utilize the mathematical foundation of Schur's decomposition theory as summarized below to ensure the existence of the solution: An arbitrary square matrix A can be transformed to a triangular matrix T as: where Q is a unitary matrix.As (4) belongs to the class of similarity transformation, matrix T has the same spectrum as matrix A .On the other hand, since T is triangular, its eigenvalues are explicitly listed as its diagonal entries 5 .This forms the mathematical basis for the use of variational quantum algorithms, which iteratively triangularize non-Hermitian matrices with unitary operators.Once a matrix has been triangularized, its eigenvalues can be readily determined by measuring the entries along the diagonal of the triangularized matrix.

Quantum process snapshot
Given the above mathematical basis, the next question to address is when one can terminate the iterative process.We develop a quantum process snapshot technique to effectively measure the proximity of an N-dimensional matrix to a triangular matrix, and it only requires a computational cost of O(log 2 N) gates.To ascertain if an N-dimensional matrix is an upper triangular one, we need to inspect its N(N − 1)/2 entries below the diagonal.On a classical computer, this is a simple task, achieved by verifying if these entries are zero.However, in quantum computation, the matrix (or operator) elements are not directly observable without measurement.Therefore, a thorough check would involve evaluating N(N − 1)/2 entries, necessitating at least O(N 2 ) operations.Such an approach would eliminate any quantum advantage we might derive from quantum computing.Inspired by the quantum process tomography technique 31 , we develop the quantum process snapshot.It captures all N(N − 1)/2 characteristics by one single quantum circuit and thus is denoted as a "snapshot" of the quantum process.
This implies the possibility of measuring |e i � from the quantum state M|e j � should be zero, for j = 1, 2, . . ., N − 1 and i = j + 1, j + 2, . . ., N .Naively, to go through all quantum states |e j � , the procedure above needs to repeat N-1 times.By leveraging quantum parallelism and quantum interference, quantum states M|e 1 �, M|e 2 �, . . ., M|e N � can be simultaneously encoded to one quantum state as shown in Fig. 1.
As this single circuit could capture all characteristics for determining its triangular property, we call it quantum process snapshot (QPS).Performing the M operator to the augmented qubits gives the quantum states ( 6): where |e i_j � = |i� ⊗ |j� , which is the (iN + j) th computational basis.
This quantum state corresponds to the vector: All matrix entry evaluations needed in (5) are then directly accessible from the measurement.For instance, the entry M 21 corresponding to the quantum measurement |�e 2 |M|e 1 �| 2 can be evaluated by checking the possibility corresponding to the vector in (8).This vector corresponds to the measurement result of quantum state |e 0_2 � as defined in (6): Therefore, the QPS method could efficiently determine whether the implemented operator M is triangular or not.However, in a quantum circuit, only the unitary operator is directly implementable, which only covers a small class of matrices (unitary matrix).This issue will be resolved in the next subsection, where we present a novel way to implement the non-unitary operator directly into the QPS process.

Non-unitary gate operation
A novel method is devised for executing non-unitary gate operations within a unitary gate-based quantum machine.This innovative approach empowers variational quantum algorithms to compute the eigenvalues of any given matrix.In order to incorporate an arbitrary non-unitary operator M into a quantum circuit, we start by decomposing it into a linear combination of unitaries, i.e., M = i c i A i , with real number coefficients c i and a set of unitary operators A i .We note that such decomposition into sum of unitaries was originally proposed by Childs and Wiebe in Hamiltonian simulations 32 .In this paper, ( 2) is assumed for such decomposition, but we note that other efficient decomposition of a matrix to unitaries do exist 33 .However, exploring other ways of decomposition is beyond the scope of this paper.
In the original VQE method, to efficiently evaluate the expectation value of a Hamiltonian by a quantum device, the matrix H is decomposed as a linear combination of the simple Pauli operators (which are also unitary).The quantum expectation evaluation is then applied to each Pauli operator term U (θ) * σ s U (θ) .However, such implementation in QPS is not enough to determine a triangular matrix.This is because an independent implementation of unitary operators could only provide the value of each (|�e l |c i A i |e m �| 2 ) .However, the summation of ( 6) As a comparison, in this subsection, we will adapt and modify the method introduced by Berry and Childs 34 , which implements the non-unitary gate operations directly.Similar ideas to implement a non-unitary operator to a quantum state |ψ� are summarized by Childs et al. 35 .As an example, consider a non-unitary operator (matrix) M = (U 1 + U 2 ) .To implement M on a state |ψ� , by adding one ancilla qubit, we start with the quantum state |0� ⊗ |ψ� .Applying the Hadamard gate to the first qubit gives 1 √ 2 (|0� + |1�) ⊗ |ψ� .We then perform the controlled unitary gate In the end, by applying Hadamard gate to the ancilla qubit, we end up with the quantum state If we measure the first qubit and obtain the |0� outcome, then we successfully prepare the qubit in the state (U 1 + U 2 )|ψ� .Therefore, although all gate operations used are unitary, the non-unitary operation can still be implemented to a quantum state with an ancilla qubit.The non-unitary operation is essentially coded under a sub-state relating to ancilla qubit.
Berry and Childs 34 and Childs and Wiebe 32 use this idea to simulate hamiltonian dynamics of quantum state ψ , where the Hamiltonian is decomposed in the form of M = 2 n −1 i=0 c i A i , where c i is a positive constant number and A i is a unitary operator.Firstly, a unitary operator B applied to the n ancilla qubits is defined as 34 : where s = 2 n −1 i=0 c i .It is also important to realize that although two possible solutions exist for each √ c i , either of them could be selected to prepare the state.They will result in the same final output which is related to , where the superscript * denotes complex conjugation.where M ′ = 2 n −1 i=0 |c i |A i and the ancillary state in state |φ� is supported in the subspace orthogonal to |0� .When the ancilla qubits are measured in state |0� , the working qubits are in state M ′ |ψ�.

Each unitary operator
From (13), it can be seen that the implemented operator M ′ only equals the attempted non-unitary operator M = 2 n −1 i=0 c i A i , when all coefficient c i is real and non-negative.However, in general, the unitary decomposition in (2) will have negative or even complex c i .To enable negative weights in the decomposition as shown (2), the controlled operator W needs to be modified as below: where the operator B c is defined as: where the superscript * denotes complex conjugate.
Apply the modified operator to the quantum state |0� ⊗ |ψ� gives: where s = 2 n −1 j=0 |c i | and the ancillary state in state |φ� is supported in the subspace orthogonal to |0� .In such implementation, instead of ( √ c i ) * • √ c i = |c i | , the coefficient corresponding to the unitary matrix A i will become This makes the implemented operator always the same as M = 2 n −1 i=0 c i A i .The mathematical proof of ( 16) can be seen in section 'Proof of non-unitary gate operation implementation' of 'Methods' .

Variational quantum universal eigensolver
We devise the Variational Quantum Universal Eigensolver (VQUE) to empower the Variational Quantum Algorithms (VQA) to evaluate eigenvalues of non-Hermitian matrices.Our approach is underpinned by Schur's decomposition theory, which assures the existence of a unitary matrix capable of triangularizing the www.nature.com/scientificreports/non-Hermitian matrix.This unitary matrix is iteratively located using the VQA.To ensure the efficiency of this search, we introduce the Quantum Process Snapshot (QPS) method.Furthermore, we have developed a novel cost function to optimize matrix triangularization.This function is efficiently evaluated using the QPS method, and it not only enhances the original cost function in the Variational Quantum Eigensolver (VQE) method but also provides an additional measure of accuracy.
Figure 2 shows the schematic diagram of the introduced variational quantum universal eigensolver.To calculate the eigenvalue of an arbitrary matrix M , firstly, the classical computer prepares the inputs to the VQUE: (1) Decompose the matrix M as a linear combination of unitaries M = i c i A i ; (2) Initialize the variational parameters θ 0 for the ansatz circuit Q(θ) .After taking the inputs, VQUE enters an optimization process to train the ansatz Q(θ) to triangularize the matrix M , i.e., where matrix A i is unitary and thus can be directly realized by quantum gate operation.
By utilizing the QPS method above, the cost function C(θ) in ( 18) is designed to quantify the triangulation extent: where each (|�e l |T|e m �| 2 ) is just the possibility of measuring the quantum state |e m_l � (see the definition in ( 6)).Thus the summation is essentially counting the number of measurements N C which fall into the quantum states set C = {|e m_l � where l = 2, 3, . . ., N and m = 1, 2, . . ., l} , compared with the total number of shots N shots , i.e.C(θ) ≡ N C /N shots .
It is also crucial to take into account that T is implemented probabilistically (see ( 16)), thus the number of successful implementations N si might be significantly less than the number of circuit executions N shots .However, insufficient N si may create statistical fluctuation to evaluate C(θ) , which makes it difficult to accurately compute its gradient.Fortunately, based on the theory proposed by Cerezo 22 , this challenge can be overcome by simply increasing a reasonable number of circuit executions, as the possibility of occurring large relative errors will exponentially reduced with a linearly increased N si .Therefore, in practice, we can measure all the qubits and repeat the experiment N shot number of times to numerically check how large N shot needs to be to achieve the required accuracy ε .In Section 'Influence of Number of Circuit Executions' of 'Parametric study' , this is dem- onstrated in the study of evaluating the eigenvalues of a 4 by 4 matrix.
From the measurement results by the quantum computer, the cost function C(θ) can be computed.The parameters θ are then adjusted by an optimizer in the classical computer (e.g., gradient descent), which is used by the ansatz in the next round of the loop.It is easy to verify that when C(θ) is approaching zero, so is (5), and thus the matrix T can be quantified as close to a triangular matrix.( 17) In theory, the target of VQUE is to find the parameters θ which minimizes the cost function C(θ) , In practice, however, the main interest is usually focused on finding a close enough solution, within a preset error threshold based on the application purpose.However, in the traditional VQE method, no quantification index is directly available, and the extent of optimization is unknown before finding the global optimum solution.This might become the main bottleneck for applying the traditional VQE technique to solve real-world problems, where the numerical error is expected to be within a specific range.Besides, even with extensive effort, phenomena like the barren plateau 36 and noise introduced in NISQ devices might still deeply influence the optimization performance, thus the optimized solution might be impossible to find.Efforts are still being paid to mitigate the barren plateau issues 37 .
This challenge is partially overcome by the cost function introduced in ( 18), whose value is designed as a quantification index of a triangular matrix.The more a matrix is triangularized, the closer its diagonal entries are to its eigenvalues.Thus a customized error threshold ε can be used as a termination condition for the opti- mization process, as follows Once the termination condition is reached with Q(θ op ) , the diagonal elements of the matrix T = Q(θ op ) † MQ(θ op ) are used as the estimate of the eigenvalues of the original matrix.Therefore, performing the Hadamard test to each unitary Q(θ op ) † A i Q(θ op ) (note M = i c i A i ) with the intial state |e i � , the real and imagery parts of i th eigenvalue can be evaluated.
It is also worth noting that the algorithm designed here is to reach out to all eigenvalues similar to the classical eigenvalue algorithm.It might be possible to adapt and modify this method if a specific subset of eigenvalues are the main interest.

Resource analysis
In this subsection, we analyze the escalation of computational resources in correlation with the expansion of the network size.Consider an N-dimensional matrix that can be decomposed into the sum of n unitaries.Quantum hardware possesses the ability to store an N-dimensional vector and perform matrix operations using only O(log 2 N) qubits, requiring exponentially fewer resources compared to classical computers.The computational complexity of executing VQUE in each iteration to triangularize the matrix approximates to O(log 2 N + n 2 ) .This complexity encompasses: 1).O(log 2 N) gates for the execution of the quantum process snapshot 2).O(n 2 ) gates for the implementation of the non-unitary gate operation 3).Additional gates required for the ansatz.Consequently, if a matrix can be represented by a limited number of unitaries (i.e., N >> n ), the proposed VQUE could potentially offer an exponential quantum advantage in matrix triangularization.However, the quantum advantage will diminish to linear to produce all N eigenvalues.It is critical to underscore that, akin to all other Variational Quantum Algorithms (VQAs), the extent of quantum advantage conferred by VQUE is largely contingent on the optimization process.

Parametric study
In this section, we showcase the numerical results from our VQUE implementation.First, to validate the devised method, we implement the VQUE on IBM's noise-free quantum simulator to assess the eigenvalues of a non-Hermitian matrix.Following this validation, we deploy the VQUE on a real IBM quantum computer (ibmq-quito, 5 qubits, 16 quantum volume).To evaluate scalability and universality, we apply the VQUE to a 1024-dimensional matrix and a non-diagonalizable matrix.To primarily emphasize the performance assessment of the algorithm, we conduct these tests on the noise-free quantum simulator.Ultimately, we employ our algorithm for the small signal analysis of an actual power network.

Verification of VQUE on by a simple illustrative example
We will demonstrate our method with a simple example via a quantum simulator and using a real quantum device.

Implementation in noise-free quantum simulator
The goal is to minimze the cost function in (18) and obtain the optimal parameters, In this subsection, we use the VQUE approach to evaluate eigenvalues of the matrix M defined as follows, Applying (2) decomposes matrix M as a linear combination of unitaries M = −2 • (σ Note there is a negative coefficient -2 in this decomposition, thus the operator M could not be directly imple- mented by the technique presented in 34 .To implement the matrix in a quantum circuit, firstly, we perform the square root operation on the coefficient vector −2 3 , giving 1.4142j −1.7321 , which is then normalized to 0.6325j −0.7746 by a constant 2.2361. To search the unitary matrix Q which triangularizes the matrix (i.e.T = QMQ † ), 5 qubits are used as shown in Fig. 3.We use the RZ-SX-RZ gates to form the variational circuit, which results in an ansatz with 4 parameters.
Firstly, the B circuit defined by (11) converts the ancilla qubit to the quantum state: In the meanwhile, the quantum augmentation circuit efficiently prepares a quantum state: which corresponds to the vectors:  and e i is a 4-dimensional unit vector with 1 in its i-th entry.
The ansatz Q(θ) is then performed to implement the matrix transformation T = Q(θ) † MQ(θ) .Subsequently, we perform the matrix formulation circuit that kickbacks the phase term in (21) to the corresponding operators (σ I ⊗ σ I ) and (σ y ⊗ σ x ) , respectively.In the end, the inverse ansatz circuit ( Q(θ) † ) and the inverse circuit of B c (i.e., B c † ) are applied to the augmented qubits and ancilla qubits, respectively, which yields the quantum state: where the ancillary state |φ� is supported in the subspace orthogonal to |0� and 5 is the coefficient s calculated from (16).Therefore, if the ancilla qubit is measured in state |0� , operator T is successfully applied to working and aug- mented qubits, which corresponds to the vector  After the cost function is evaluated, "Constrained Optimization By Linear Approximation optimizer" 38 is used to update the parameters θ , which will then be used in the ansatz Q(θ) in the next iteration.Figure 3 plots the cost function versus the number of iterations.It shows that the cost function is consistent below 0.01 after around 140 iterations.Thus the θ after 155 iteration is selected as θ op .
The Hadmard test is then applied with the optimized ansatz Q(θ op ) to evaluate the i th eigenvalue.The evaluated eigenvalues are shown in the results part in Fig. 3. Compared with the theoretical eigenvalues, 3 has the largest relative error of 0.39% , which is quite close to the theoretical one.Therefore, we conclude that: • In theory, the devised method has the capability to evaluate the eigenvalues of a general matrix with complex eigenvalues.• As the absolute error is similar for all eigenvalues, the evaluation of smaller eigenvalues might have larger relative errors.

Influence of number of circuit executions
As discussed in section 'Variational quantum universal eigensolver' of 'results' , the number of circuit executions N shots plays a key role in accurately evaluating the eigenvalues.As pointed out by Cerezo 22 , the possibility of resulting large evaluation errors should exponentially decrease with a linear increased N shots .To confirm this theoretical analysis, we make a numerical study in the same quantum circuit as shown in Fig. 3.
In the experiment, we set the accuracy threshold as the maximum error in the eigenvalue evaluations is less than 5% .From the experiment, the algorithm starts to make meaningful estimation with N shots = 300 .Thus we increased the number of shots N shots from 300 to 2800, with a step of 25.For each N shots , we repeat 100 times and count the number of results which gives the eigenvalues inside the accuracy threshold.The experiment result is shown in Fig. 4.
It can be observed that with the linearly increased number of N shots , the possibility for VQUE to make estima- tions with maximum relative error larger than 5% is close to exponentially decreased.Therefore, it is reasonable to expect overcoming the challenge of T is implemented probabilistically by an acceptable increase of N shots .

Implementation on real quantum machine
In this subsection, we implement example 1 in a real quantum machine to study its performance in the current Noisy Intermediate-Scale Quantum (NISQ) environment.We use the IBM quantum computer ibmq_manila, which is a 32-quantum volume machine with 5 qubits.Its structure is shown in the quantum machine configuration part in Fig. 5.
As shown in Fig. 5, due to the limited connectivity of the real quantum computer, the quantum entanglement between the ancilla qubit and augmented qubits will require a significant number of CNOT gates to realize in the quantum computer.Therefore, we made an engineering trick here by removing the ancilla qubit and only performing σ y ⊗ σ x to the augmented qubits as shown in the VQUE circuit in theory in Fig. 5.Such treatment is based on the fact that Q † IQ = I for an arbitrary unitary matrix Q .Therefore, searching for the ansatz to (24) The circuit is then compiled to fit the configuration of the quantum machine as shown in Fig. 5.It can be observed that the noise introduced by the quantum computer will significantly influence the classical optimizer.As a result, the cost function can only converge to a local minimum of around 2, which results in the estimated eigenvalues having a relative error from 14.38% to 25.42%.

VQUE scalability demontration
In this subsection, the VQUE algorithm estimates the eigenvalues of a 1024-dimensional matrix in the noise-free quantum simulator.The matrix is composed of two unitary matrices:   www.nature.com/scientificreports/

Verification of VQUE for non-diagonalizable matrix
In this subsection, to demonstrate the universality, the VQUE is applied to evaluate the eigenvalue of a nondiagonalizable matrix N: It is easy to verify that this matrix only has three linear independent eigenvectors and thus is not similar to any diagonal matrices 5 .Decomposing the matrix by (2) gives the 13 coefficients, which is significantly more complicated than the matrix used in the previous sections.As shown in Fig. 7, when the same ansatz as in previous cases is used, the cost function could only converge around 4, which is too far away from a triangular matrix.By this ansatz, the evaluated eigenvalues are shown on the left side of the results in Fig. 7.It can be observed that 3 has the largest evaluation error of ǫ = 1.3079 − 0.3270 , which is 134.82% of its original eigenvalue.This indicates that the original eigenvalue 3 is submerged by the error, which is consistent with the conclusion we have drawn from the cost function.
To precisely evaluate the eigenvalues, we slightly modify the ansatz as shown in the "C-Ansatz" of Fig. 7.This complicated ansatz enables the VQUE to search for the solution in a much larger subspace.With the empowered ansatz, the cost function could then be brought down to 0.001 after 150 iterations.By this optimized ansatz, the eigenvalues are listed on the right side of the results in Fig 7 .At this time, the largest relative error appears in eigenalue 1 , which is 7.37% of original eigenvalue.Compared the results with the one from the simple anstanz in the beginning, the eigenvalue evaluation accuracy is dramatically improved from 134.82% error to 7.37% error.
It is essential to realize that it is our devised method that enables such a trial-and-error process.In the original version of VQE, only the smallest eigenvalue is available, and there is no indication of how close the evaluated eigenvalue is to the theoretical one.This is critical in practical applications, where the numerical error is expected to be controlled inside a specific range.However, in today's NISQ machines, the variational method is not necessarily guaranteed to work.For instance, the barren plateau phenomenon will require exponential precision to train the ansatz.Besides, a universal parameterized circuit itself requires a very complex structure with a tremendousnumer of parameters (e.g. 39introduced a depth 7, 24 parameters universal ansatz for 2 qubits).Such complexity will not only be deeply influenced by the implementation noise but also increase the possibility of converging to a sub-optimum during non-convex optimization 23 .
In the standard VQE approach, regardless of the ansatz chosen, the algorithm consistently yields an eigenvalue estimate.However, it does not provide insights into the accuracy of the predicted eigenvalue, necessitating the use of a sufficiently comprehensive ansatz to guarantee a dependable evaluation.In contrast, the cost function within our proposed methodology serves a dual purpose, also acting as an indicator of accuracy.As shown in example 3, this feature allows for the exploration of a shallower tentative ansatz through a process of trial and error, as it offers immediate feedback on the reliability of the results.Our introduced method thus opens an engineering window to solve the problem without explicitly overcoming all the challenges above.A threshold value can be set based on the specific application.We could then start with a simple ansatz and gradually increase its complexity to attempt to reduce the C(θ) below the threshold value.If a solution that matches the requirement Q.E.D.
then implemented by the controlled operator select(A) defined as: Therefore, applying the unitary operator W = (B † ⊗ I)(select(A))(B ⊗ I) to the state |0� ⊗ |ψ� gives 34 :
σ I ⊗ σ I ⊗ σ I ⊗ σ I ⊗ σ I ⊗ σ I ⊗ σ I ⊗ σ I ⊗ σ I ⊗ σ z and σ I ⊗ σ I ⊗ σ z ⊗ σ I ⊗ σ x ⊗ σ I ⊗ σ I ⊗ σ z ⊗ σ I ⊗ σ I , with coefficients −1 and 2 respectively.An entangled RX-RY gates ansatz as shown in Fig. 6 is selected as the parameterized circuit.As shown in the training process, the cost function is converged to 20, which implies the average value below the diagonal elements is around 3.38184 • 10 −5 .It has 1024 eigenvalues in the values set of {−1 + 2j, −1 − 2j, 1 + 2j, 1 − 2j} .The estimated eigenvalues are shown in the results part of Fig. 6.It can be observed that all the estimated eigenvalues are within the relative error of 2% , which maintains a very high estimation accuracy.This indicates that VQUE has outstanding scalability and thus brings a huge potential to resolving ultra-large size eigenvalue problems.